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ABSTRACT 

Algorithms  for  efficient  generation  of  random  numbers 
from  various  probability  distributions  are  presented,  in  both  a 
flowchart  form  and  as  a  sample  Fortran  subroutine.  Twenty- 
two  different  distributions,  including  all  commonly  encountered 
discrete  and  continuous  functions,  the  Weibull,  Johnson,  and 
Pearson  families  of  empirical  distributions,  and  histogram  dis¬ 
tributions,  are  covered.  The  general  techniques  to  apply  in 
deriving  a  random  number  selection  scheme  for  an  arbitrary 
distribution  are  discussed.  A  machine- independent  subroutine 
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EXECUTIVE  SUMMARY 


Monte  Carlo  simulation  is  one  of  the  most  powerful  and  commonly 
used  techniques  for  analyzing  complex  physical  problems.  Applications  can 
be  found  in  many  diverse  areas  from  radiation  transport  to  river  basin 
modeling.  Important  Navy  applications  include  analysis  of  antisubmarine 
warfare  exercises  and  operations,  prediction  of  aircraft  or  sensor  perform¬ 
ance,  tactical  analysis,  and  matrix  game  solutions  where  random  processes 
are  considered  to  be  of  particular  importance.  The  range  of  applications  has 
been  broadening  and  the  size,  complexity,  and  computational  effort  required 
have  been  increasing.  However,  such  developments  are  expected  and  de¬ 
sirable  since  increased  realism  is  concomitant  with  more  complex  and  exten¬ 
sive  problem  descriptions. 

In  recognition  of  such  trends,  the  requirements  for  improved  simula¬ 
tion  techniques  are  becoming  more  pressing.  Unfortunately,  methods  for 
achieving  greater  efficiency  are  frequently  overlooked  in  developing  simula¬ 
tions.  This  can  generally  be  attributed  to  one  or  more  of  the  following  reasons: 

•  Analysts  usually  seek  advanced  computer  systems  to  perform 
more  complex  simulation  studies  by  exploiting  increased 
speed  and/or  storage  capabilities.  This  is  often  achieved 

at  a  considerably  increased  expense . 

•  Many  efficient  simulation  methods  have  evolved  for  specialized 
applications.  For  example,  some  of  the  most  impressive 
Monte  Carlo  techniques  have  been  developed  in  radiation  trans¬ 
port,  a  discipline  that  does  not  overlap  into  areas  where  even 

a  small  number  of  simulation  analysts  are  working. 

•  Known  techniques  are  not  developed  to  the  point  where  they  can 
be  easily  understood  or  applied  by  even  a  small  fraction  of  the 
analysts  who  are  performing  simulation  studies  or  developing 
simulation  models. 
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1.  INTRODUCTION 


In  developing  any  Monte  Carlo  simulation,  it  is  necessary  to  generate 
random  numbers  from  the  stochastic  models  used.  In  Volume  I,  the  process 
and  techniques  of  selecting  probability  models  for  the  simulation  were  pre¬ 
sented.  The  objective  of  this  volume  is  to  provide  a  convenierV  source  of 
efficient  and  simple  random  number  generators  for  all  the  probability  dis¬ 
tributions  considered  in  Volume  I.  To  this  end  flow  charts  and  FORTRAN 
listings  of  these  random  number  generators  are  provided  here  as  well  as 
descriptions  of  the  techniques  employed. 

It  is  the  purpose  of  this  document  to  provide  a  convenient  mechanism 
to  select  and  implement  these  random  number  generators  without  having  to 
resort  to  an  understanding  of  the  underlying  concepts  used  in  their  develop¬ 
ment.  Accordingly,  the  remainder  of  this  report  has  been  organized  as 
follows: 

•  SECTION  2,  "Efficiency  Comparison  of  Random  Number 
Generators,"  demonstrates  improvements  in  running  times 
expected  from  using  the  techniques  developed  here  over  those 
commonly  used.  This  section  has  been  included  to  provide  an 
appreciation  for  the  magnitude  of  improvements  possible  in 
using  the  techniques  described  herein. 

•  SECTION  3,  "Generation  of  Random  Numbers  from  Selected  Dis¬ 
tributions,"  provides  algorithms  defined  by  flow  diagrams  and 
standard  Fortran  subroutines  that  can  be  applied  directly.  This 
section  is  introduced  with  a  convenient  summary  table  defining 
where  in  the  section  a  specific  algorithm  can  be  found. 

•  appendix  A,  "Fundamental  Considerations  for  Generation  of 
Random  Numbers,"  describes  the  fundamentals  on  which  random 
number  generation  techniques  for  arbitrary  distributions  can  be 
developed. 
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•  Appendix  B,  "MIRAN  -  A  Machine  Independent  Package  For 
Generating  Uniform  Random  Numbers,  ”  describes  a  uniform 
random  number  generator  that  can  be  used  on  any  machine 
that  does  not  have  a  reliable  generator  or  on  several  different 
machines  where  identical  random  numbers  are  to  be  generated 
for  comparison  and  cross  checking. 

Before  proceeding  it  must  be  recognized  that  a  "good”  uniform  ran¬ 
dom  number  generator  is  generally  assumed  to  be  available  to  the  user. 

This  is  often  not  the  case,  although  most  computerc  today  have  uniform 
random  number  generators  included  as  part  of  the  system  software.  Un¬ 
fortunately,  many  of  the  uniform  random  number  generators  in  current 
use  do  not  adequately  approximate  randomness  to  be  sufficient  for  all  Monte 
Carlo  calculations.  To  alleviate  this  difficulty,  a  machine  independent 
package  for  generating  uniform  random  numbers  is  provided  (Appendix  B). 
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2.  COMPARISON  OF  RANDOM  NUMBER  GENERATION  PROCEDURES 


The  improvements  in  calculational  efficiency  realized  by  using  the 
random  number  generation  techniques  provided  here  depend  on  the  particular 
problem.  However,  by  utilizing  these  techniques,  near  optimum  results  can 
be  assured. 

It  is  of  interest  to  compare  the  random  number  generation  techniques 
presented  here  with  those  commonly  used  to  generate  random  numbers.  This 
comparison  was  performed  during  the  course  of  the  study  for  several  distri¬ 
butions,  and  it  was  found  that  improvements  in  computer  time  of  factors  vary¬ 
ing  from  2  to  5  were  possible.  Results  for  a  few  of  the  more  common  distri¬ 
butions  are  shown  in  Table  2. 1  which  compares  the  running  times  of  the 
preferred  techniques  with  those  commonly  used.  For  example,  consider  the 
normal  (or  Gaussian)  distribution.  The  usual  procedure  is  to  generate  12 

random  numbers  uniformly  distributed  over  the  interval  [0,  ll  say  Rj . R^. 

and  determine 

12 

Rn  ■  Z  V6  • 

i=l 

/c\ 

By  virtue  of  the  central  limit  theorem, '  '  R^  is  approximately  distributed 

according  to  the  normal  distribution.  Assembly  language  time  on  a  Univac 

1108  was  105  microseconds  per  calculation  using  this  approach.  Procedures 

studied  here  were  the  rejection  technique  (see  Appendix  A)  and  a  technique 

(5) 

developed  by  Marsaglia.  The  corresponding  running  times  were  respec¬ 
tively  74  and  30  microseconds.  Not  only  are  the  running  times  significantly 
reduced,  but  also  the  more  efficient  ones  presented  here  are  exact  (within 

machine  roundoff  errors). 

Similar  results  were  obtained  with  the  exponential  distribution  where 
the  Marsaglia  technique  gave  a  reduction  in  running  times  of  a  factor  greater 
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than  three  (Table  2.1).  The  standard  method  used  is  the  inverse  (see 
Appendix  A).  The  rejection  method  is  discussed  in  Appendix  A  and  the 
Marsaglia  method  is  reported  in  Ref.  3. 

As  implied  above,  there  are  several  methods  that  may  be  used  to 
generate  random  numbers  for  a  given  distribution.  However,  where  alternate 
approaches  could  be  identified  or  developed,  comparisons  were  made  and  the 
most  efficient  procedure  selected.  These  generators  are  presented  in  the 
next  section. 

It  should  be  noted  that  the  more  efficient  techniques  are  slightly  more 
complex  to  program;  however,  the  slight  additional  effort  involved  gener¬ 
ally  pays  off  substantially  in  computer  time. 


TABLE  2  .1 

Running  Time  Comparisons  Random  Number  Generators  For 
The  Normal  and  Exponential  Distributions11 


Commonly 

Marsagliaa 

Used 

Rejection 

Distribution 

Technique 

Technique 

Technique 

Exponential 

64 

29 

19 

Normal 

(Gaussian) 

105 

74 

30 

aSee  Appendix  A  for  a  brief  description  of  these  techniques, 
b 

All  times  in  microseconds  of  UNIV^C  1108  Assembly  Language  time. 
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3.  GENERATION  OF  RANDOM  NUMBERS  FROM  SELECTED 

DISTRIBUTIONS 


In  this  section,  efficient  algorithms  are  presented  for  a  large  number 
of  probability  distributions.  These  are  summarized  in  Table  3-1  which 
gives  the  name  of  the  distribution,  the  theoretical  form,  parameters  in  the 
distribution  to  be  specified  by  the  user,  other  random  number  generators 
used,  and  where  the  particular  routines  or  algorithms  can  be  found  in  this 
section  of  the  report.  Also  shown  under  the  name  of  the  distribution  is  the 
FORTRAN  subroutine  name  assigned  to  the  random  variable. 

Once  a  distribution  of  interest  has  been  identified,  it  is  only  necessary 
to  define  the  values  of  the  parameters  indicated  and  to  implement  the 
algorithm  from  the  specified  pages  of  this  section.  In  the  subroutines, 
the  parameters  are  represented  by  mnemonics  which  should  be  recog¬ 
nizable.  For  example,  SIG  is  used  to  represent  a  and  SIGSQ  to  repre- 

2 

sent  c r  .  in  some  places  the  mnemonic  starts  with  an  A  to  provide  a  float¬ 
ing  point  value  such  as  ALAM  for  X  . 

It  will  be  noted  that  certain  distributions  rely  on  other  distributions 
to  generate  random  numbers.  For  example,  generation  of  random  numbers 
for  the  Rayleigh  distribution  requires  random  numbers  from  an  exponential 
distribution.  The  exponential  distribution  in  turn  depends  on  a  uniform 
random  number  generator.  Based  on  the  frequent  requirement  for  the  uni¬ 
form,  exponential  and  normal  distribution,  it  is  usually  convenient  to  pro¬ 
vide  a  basic  random  number  generation  package  consisting  of  subroutines 
to  generate  uniform,  exponential,  and  normal  random  variables  as  an  inte¬ 
gral  part  of  any  complex  simulation  program.  Throughout  this  section  these 
three  random  number  generation  subroutines  will  appear  as  UNFRN(R), 
EXPRN(R),  and  ANRMRN(R),  respectively,  where  R  is  a  dummy  function 


TABLE  3.1 


Efficient  Algorithms  for  a  Large  Number  of  Probability  Distributions 


Namr  of 
Distribution 
(Function  Title) 

Functional  Form 

Parameter* 

To  Be  Specified 

Other 

Random  Number 
Generator*  U**d 

Location  of 

Algorithm  to  Generate 
Random  Number* 

Subjection 

Page 

Uniform 

(UNFRN) 

b^;  *sxsb 

a,  b 

None 

3.1 

10 

Exponential 

(FXPRN) 

1,A  e-K«’*)A]  ;  x>« 

A  >0 

A,f 

Uniform 

3.2 

12 

Normal 

(ANMRN) 

1  -(x-„)2/2o2 

ay[r. 

P>  e 

Uniform, 

Exponential 

3.3 

14 

Binomial 

(KWNOM) 

(kn)pk(l-p  >n'k; 

k  «  0,1,  ....  n 

n»  P 

Uniform, 

Exponential 

3.4 

17 

Multinomial 

(MULNOM) 

(k1kJ...km)pi1pJV--Pmm 

pl  4  •  *  •  *  pm  *  1 
k1«kJ4...+k|n.,» 

"•  p| . Pn 

Uniform 

3.5 

22 

Poisson 

(KPOIS) 

-A  Ak 

e  irr;  *>0 

k  *  0, 1, . . , 

A 

Uniform 

3.6 

24 

Hyper -geometric 
(KHYPRO) 

(?)  «■).„„ 
m  ■ 

k  *  0,  1 . M 

M,  N,  n 

Uniform 

3.7 

26 

Geometric 

(KGEOM) 

puV'1 

k  •  1,2,3,... 

P 

Uniform, 

Exponential 

3.8 

■ 

Pascal  (also 
called  negative 
binomial) 

(KPASCL1 

(n  +  k'1)  u-P)nPk 

k  >0,1 . n 

n,P 

Uniform, 

'X.7pcn*ntial 

3.9 

31 
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TABLE  3.1  (Continued) 
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Distribution 

Parameter* 

O*iior 

Random  Number 

Location  of 
Algorithm  to  Generate 
Random  Numbers 

j  (Function  Title) 

Functional  Form 

To  Be  Specified 

Cenentors  Uxcd 

Subsect  on 

Tape 

We  i  bull 
(V.IBLR.V) 

itx-O’ 

.isiill 
-1.  * 

’ll  * 

Expo:  ntul 

3.19 

53 

*  2  < 

p?,A  >  0 

1 

n 

3.  ?C  1 

55 

L  N2- 

(x-<) 

•  •  7 .  y 

(SLRX) 

»*? 

x>  0 

X  *  • 

SB:  ^ 

A 

(X  -  *  >  A  -  >  *  0 

i,  y. 

Normal 

3.20.2 

57 

exp 

VI 

|  2  1*'*  1  IM  A-x.i  j 
1  !),X  >0 

(  s  x  s  *  ♦  A 

J  ) 

V  -7= 
u  \  2* 

1 

1  (x  ♦  <)J  - 

1,  >■.*.' 

Normal 

3.20.3 

59 

:sv?.:r, 

np 

'2'  »‘f"|(5r) 

2  x1/2!VJ 

))] 

i),X  >0 

Pearson 

System 

Type  I- 

m,  m2 

Carr.-ra 

3.21.1 

■1 

(TYP’.RNl 

C(" 

■ 

Type  R: 

a,  m 

Gamrr.x 

3.21.2 

03 

{TYT2R.S1 

C(l* 

2  \m 

xM  m  >  -1 

»2/  -.«*«» 

Type  nt: 

3.21.3 

65 

(TYP2RN") 

c(l. 

_\  y»  -rx  -x  <  x  <  x 

l) 

F.  a 

Gamma 
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Kau*e  td 

Diet  nl««t  ion 

Functional  form 

Parameter* 

To  ke  (pectfled 

(ItVr 

Random  Number 

Generators  UMd 

Loco 
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argument.  In  the  flow  diagrams,  these  are  indicated  as  U(0, 1),  E(0, 1) 
and  N(0, 1),  respectively. 

3. 1  UNIFORM  RANDOM  NUMBER  GENERATORS 

The  uniform  random  number  generator  is,  of  course,  fundamental 
to  all  random  number  generation.  For  the  purposes  here,  it  is  assumed 
that  the  computer  system  available  will  have  such  a  generator  as  part  of 
the  basic  software  package.  If  one  is  not  available  or  the  generator  is 
expected  to  be  faulty,  the  machine  independent  package  presented  in  Ap¬ 
pendix  B  (MIRAN)  can  be  used.  The  following  paragraphs  describe  the 
technique  used  in  most  computers  for  generating  random  numbers  and  pro¬ 
vide  insight  into  the  assessment  of  such  generators. 

The  method  used  for  almost  all  uniform  generators  is  the  multiplica- 

(7) 

tive  congruential  method. v  A  sequence  of  integers,  x  ,Xj, . . . ,  is  generated 
by  the  congruence 

p 

xn+l  =  xn'*(mod2  )  . 

Here  P  is  the  number  of  bits  (excluding  sign)  in  a  word  on  the  particular 
computer  employed  and  X  is  called  the  generator  which  is  a  carefully  selected 
integer  as  described  below.  From  this  sequence  random  fractions  are  pro¬ 
duced  using 

R  =  x  .2‘P  . 
n  n 

The  sequence  of  random  fractions,  R^Rg, . . . ,  is  output  by  the  subroutine  in 
floating  point  form. 

On  most  computers  the  multiplicative  congruential  method  is  accom¬ 
plished  by  an  integer  multiplication  of  x  and  X.  Only  the  low-order  half 

n 

(P  bits)  of  the  product  is  retained  as  xn+j*  This  is  then  treated  as  a  binary 
fraction,  converted  to  floating  point,  and  normalized. 
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This  method  is  fast  and  will  produce  numbers  whose  properties  ap¬ 
proximate  randomness  sufficiently  close  for  valid  use  in  Monte  Carlo 
simulations  provided  the  following  caveats  are  observed. 

1.  Choose  a  generator,  X  ,  with  particular  care.  In  particular, 

generators  with  a  small  number  of  T  bits  in  their  binary  repre¬ 
sentation  should  be  avoided.  A  number  of  generators  of  the  form 
2^  ±  3,  2^4  ±  3,  §16  ±  ^  etc. ,  are  particularly  abundant.  At 
one  time,  they  were  used  because  they  were  thought  to  be  good 
and  especially  fast.  However,  further  research  has  shown  them 
to  be  faulty  and  a  number  of  simulations  have  produced  erroneous 
results  as  a  consequence.  Small  generators  such  as  X  =  101  p 

are  also  faulty  and  must  be  avoided.  The  generators  X  *  515  or  x  =  5 
have  been  well  tested  and  are  quite  safe  to  use. ^ 

2.  Check  the  computer  word  length.  It  is  best  for  P  to  be  at 
least  35  in  the  congruence.  For  machines  with  P  <  32  a  multi¬ 
ple  precision  multiplication  should  be  used  to  generate  an  ade¬ 
quate  congruence. 

3.  Do  not  trust,  on  blind  faith,  random  number  routines  distributed 
by  the  computer  manufacturers  with  standard  subroutine  libraries. 
These  have  been  found  to  contain,  with  high  probability,  the  faulty 
generator  values. 

The  uniform  random  number  generator  will  be  referred  to  as  UNFRNfR^ 
in  subsequent  routines  and  U(0, 1)  in  the  flow  diagrams. 
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3. 2  EXPONENTIAL  DISTRIBUTION 


The  simplest  method  to  generate  random  numbers  from  the  exponential 
distribution,  f(x)  =  e  ,  is  to  use  the  inverse  solution, 

x  =  -lnfru)  , 

where  R .  is  a  uniform  random  number.  This  is  not,  however,  the  fastest 

U  (3) 

method.  An  extremely  rapid  technique  has  been  de^  'loped  by  G.  Marsaglia 

which,  although  it  is  several  times  faster  than  the  logarithm,  requires  a 

sizable  block  of  computer  storage  (—600  words).  When  computer  storage 

is  critical  or  when  the  exponential  distrib  ion  is  not  of  crucial  importance, 

the  Von  Neumann  rejection  technique  is  a  good  general  method.  This  method, 

usually  faster  than  the  logarithm,  is  shown  in  Fig.  3-1. 

To  select  from  a  generalized  exponential,  (l/X)e”^x~c^^,  it  is 
merely  necessary  to  select  from  e"x  then  multiply  by  X  and  add  c .  For 
best  efficiency  in  general,  the  basic  exponential  subroutine  should  select  from 

-X 

e  ,  and  it  should  be  left  up  to  the  calling  program  to  supply  the  multiplication 
and  addition  where  needed. 

The  exponential  distribution  is  referred  to  as  EXPRN(R)  in  subsequent 
routines  and  as  E(0, 1)  in  the  flow  diagrams. 

Sample  Routines 

Simplest  method  (use  inline  in  calling  program) : 

R  =  -ALOG  (UNFRN(R)) 

Von  Neumann  rejection  technique: 


FUNCTION  EXPRN(DUMMY) 
I  =  0 

100 

X  =  UNFRN(X) 

105 

Y  =  UNFRN(X) 

IF  (X.  LT.  Y)  GO  TO  120 

110 

X  =  UNFRN(X) 

IF  (X.LT.Y)  GO  TO  105 

115 

I  =  1+1 

GO  TO  100 

120 

EXPRN  =  X+I 

RETURN 

END 
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Generate  y  «-U(0, 1) 


Figure  3-1.  Random  number  generation  algorithm 
for  exponential  distribution 


3. 3  NORMAL  DISTRIBUTION 


The  normal  distribution,  f(x)  =  l/(a«/2ir)e“'x"^  '  ,  has  received 

considerable  attention  by  the  designers  of  random  number  generators.  One 

of  the  earliest  methods,  which  is  still  found  frequently  in  simulations  today, 

uses  the  central  limit  theorem  to  approximate  the  normal  by  summing  up 

(6) 

several  uniform  random  variables.  This  approach  has  two  serious  defects. 

First,  it  is  only  an  approximation.  Second,  it  is  much  slower  than  other 

methods.  The  fastest  method  by  far  is  a  technique  designed  by  G.  Marsaglia. 

However,  considerable  storage  is  needed  for  this  technique.  Another 

(4) 

technique  by  Marsaglia,'  '  illustrated  in  Fig.  3-2,  is  fairly  fast  without 
requiring  much  computer  storage.  Thia  is  the  best  technique  known  for 
general  usage. 


(5) 


As  with  the  exponential  routine,  the  basic  normal  random  number 

generator  should  be  written  to  select  from  the  normal  distribution  with  unit 

mean  and  zero  standard  deviation  (referred  to  as  ANRMRN  in  the  routines 

and  as  N(0, 1)  in  the  flow  diagrams).  It  is  then  left  up  to  the  calling  program 

to  multiply  by  the  standard  deviation  and  add  the  mean  if  a  generalized  normal 

2 

deviate  is  required.  That  is,  for  a  distribution  with  mean  4  and  variance  a  , 

the  correct  random  number  would  be  aN(0, 1)  +  p,  where  N(0, ).)  is  a  ran- 

2 

dom  number  from  a  distribution  with  4  =  0  and  cr  =  1 . 

Sample  Routine 

FUNCTION  ANRMRN  (DUMMY) 

R  =  UNFRN(R) 

IF  (R  GT.  0.  8638)  GO  TO  10 

ANRMRN  =  2.  *(UNFRN(X)  +  UNFRN(Y)  +  UNFRN(Z)  -1.5) 

RETURN 

10  IF  (R  GT.  0. 9745)  GO  TO  20 

ANRMRN  =  1. 5*(UNFRN(X)  +  UNFRN(Y)  -  1.0) 

RETURN 
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20  IF  (R.  GT.  0.  997302039)  GO  TO  100 
25  X  =  6.  *UNFRN(X)  -  3. 0 
Y  =  0.  358*UNFRN(X) 

XSQ  =  X*X 

GX  =  17.  49731196*EXP(-XSQ*.  5) 

AX  =  ABS(X) 

IF  (AX.  GT.  1. 0)  GO  TO  30 

IF  (Y.  GT.  (GX-17.  44392294  +  4. 73570326*XSQ  +  2. 15787544* AX)) 
GOTO  25 
ANRMRN  =  X 
RETURN 

30  AX3  =  2.  36785163*(3-AX)**2 

IF  (AX.  GT.  1.  5)  GO  TO  40 

IF  (Y.  GT.  (GX-AX3-2. 15787 544*(1.  5- AX)))  GO  TO  25 

ANRMRN  =  X 

RETURN 

40  IF  (Y.  GT.  (GX-AX3))  GO  TO  25 
ANRMRN  =  X 
RETURN 

100  X  =  SQRT  (9+2*EXPRN(X)) 

IF  (UNFRN(X).  GT.  3/X)  GO  TO  100 

IF  (UNFRN(X).  GT.  0.  EJ  X  =  -X 

ANRMRN  =  X 

RETURN 

END 
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3. 4  THE  BINOMIAL  DISTRIBUTION 

The  binomial  distribution,  =  (£)pk  (l-p)n-k,  is  a  discrete  distri¬ 
bution  describing  the  number  of  successes  encountered  in  a  series  of  Bernoulli 
trials.  It  has  two  parameters,  p,  the  probability  of  success  in  a  single  trial, 
and  n,  the  number  of  trials  in  the  series. 

The  algorithm  for  selection  from  the  binomial  distribution  is  divided 
into  three  subranges  for  the  parameter  p.  For  moderate  values  of  p,  the  ran¬ 
dom  number  generation  is  based  on  a  straightforward  simulation  of  the  under¬ 
lying  basis  for  the  distribution;  n  Bernoulli  trials  are  generated  and  the  num¬ 
ber  of  successes  are  counted.  For  small  values  of  p,  it  becomes  more  efficient 
to  use  a  technique  based  on  the  geometric  distribution.  Conversely,  for  large 
values  of  p  it  is  efficient  to  reverse  the  geometric  technique  and  perform  the 
counting  on  the  number  of  failures  rather  than  successes. 

For  large  values  of  n ,  all  three  algorithms  become  inefficient; 
the  computing  time  involved  is  directly  proportional  to  n  .  The  binomial 
distribution  approximates  a  normal  distribution  with  mean  np  and 
standard  deviation  ynp(l-p)  for  large  n  .  One  should  consider  replacing 
the  binomial  with  the  approximate  normal  for  large  values  of  n  (n  >  10  p/(l-p) 
or  n  >  10  (l-p)/p). 

Sample  Subroutines 

For  p  <  . 25 

FUNCTION  KBINOM  (N,  ALNQ) 

C  ALNQ  IS  -A LOG  (1.  -P) 

KBINOM  =  0 
M  =  0 

5  R  =  EXPRN(R) 

J  =  1  +  R/ALNQ 
M  -  M  +  J 

IF  (M  -  N)  10,  15,  20 
10  KBINOM  =  KBINOM  +  1 

GO  TO  5 

15  KBINOM  =  KBINOM  +  1 

20  RETURN 

END 
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For  .  25  $  p  .  75 

FUNCTION  KBINOM  (N,  P) 

KBINOM  =  0 
DO  15  M  =  1,  N 
R  =  UNFRN  (ft) 

IF  (R.  LT.  P)  KBINOM  =  KBINOM  +  1 
15  CONTINUE 
RETURN 
END 

For  p  >  .  75 

FUNCTION  KBINOM  (N,  ALNP) 

C  ALNP  IS  -A LOG  (P) 

KBINOM  =  N 
M  =  0 

5  R  =  EXPRN  (R) 

J  =  1  +R/ALNP 

M  =  M  +  J 
IF  (M-N)IO,  15,  20 
10  KBINOM  =  KBINOM  -  1 
GO  TO  5 

15  KBINOM  =  KBINOM  -1 
20  RETURN 
END 
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pk  (1  -  p)n"k 


« 


For  0. 25  :S  p  £  0. 75 


Figure  3-4.  Random  number  generation  algorithm 
for  binomial  diatribution  (continued) 


Figure  3-5.  Random  number  generation  algorithm 
for  binomial  distribution  (continued) 


3.5  THE  MULTINOMIAL  DISTRIBUTION 
The  multinomial  distribution, 


(  B  \kl  *2  k. 

p<kl>  k2’"”km)  "  ,"km/pl  p2 . pn 


is  a  generalization  of  the  binomial  distribution  to  trials  having  m  different 
outcomes  with  discrete  probabilities.  Random  number  generation  is  accom¬ 
plished  by  a  straightforward  simulation  of  the  underlying  process  of  identical 
trials.  Note  that  a  'random  number'  for  this  distribution  is  an  array  con¬ 
taining  the  number  of  realizations  of  each  possible  outcome. 

Sample  Routine 

SUBROUTINE  MULNOM  (N,  M,  K,  P) 

DIMENSION  K  (M),  P(M) 

C  P  IS  INPUT  ARRAY  OF  PROBABILITIES 
C  K  IS  OUTPUT  ARRAY  OF  OUTCOMES 
DO  10  J  -  1,  M 

10  K(J)  *  0 

DO  30  I  =  1,  N 
R  =  UNFRN  (R) 

DO  20  J  =  1,  M 

?F  (R.  LT^O )  GO  TO  30 
20  CONTINUE 

30  K(J)=K(J)  +  1 

RETURN 
END 
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Figure  3-6.  Random  number  generation  algorithm  for  multinomial 
distribution 
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3.6  POISSON  DISTRIBUTION 

-A  k 

The  Poisson  distribution,  =  e"  ,  is  a  discrete  distribution 
describing  die  number  of  occurrences  in  an  interval  when  the  rate  of  occur¬ 
rence  is  a  constant.  The  technique  for  selecting  from  the  Poisson  distribu¬ 
tion  is  a  combination-transformation  method  described  in  Ref.  2. 

The  computer  time  spent  in  tills  selection  is  directly  proportional  to 
A,  the  mean  value  of  the  Poisson  variable.  For  large  A,  this  selection 
can  be  very  time  consuming.  It  is  possible  to  approximate  the  Poisson  dis¬ 
tribution  by  %  normal  distribution  with  a  mean  of  A  and  a  standard  deviation 
of  /A  for  A  sufficiently  large  (  A  >10). 

Sample  Routine 

FUNCTION  KPOIS  (EXP LAM) 

C  EXPLAM  IS  EXP  (-LAMBDA) 

Y  =  1.0 
KPOIS  =  0 

5  Y  =  Y  *  UNFRN  (Y) 

IF  (Y.  GT.  EXPLAM)  GO  TO  10 
KPOIS  =  KPOIS  +  1 
GO  TO  5 

10  RETURN 

END 
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Figure  S-7.  Random  tiu 
Poisson  die 


Generate  R  «-U(0, 1) 


generation  algorithm  for 


3.7  HYPERGEOMETRIC  DISTRIBUTION 


The  hypergeometric  distribution, 


describes  sampling  without  replacement.  It  has  die  parameters  N  ,  the 
size  of  die  total  population,  n ,  the  size  of  the  population  sampled,  and  M , 
die  number  of  events  in  the  total  population.  The  random  variable  k  is 
the  number  of  events  occurring  in  the  sample.  The  hypergeometric  dis¬ 
tribution  is  generated  by  simulating  sampling  without  replacement. 

Sample  Routine 

FUNCTION  KHYPRG  (NTOT,  MTOT,  N) 

C  NTOT  IS  TOTAL  POPULATION  SIZE,  MTOT  IS  TOTAL 
C  EVENTS  IN  POPULATION,  N  IS  SAMPLE  SIZE 
KHYPRG  =  0 
EM  =  MTOT 
EN  =  NTOT 
DO  10  I  =  1,  N 
P  =  EM/EN 
R  =  UNFRN  (R) 

IF  (R.  GT.  P)  GO  TO  10 
KHYPRG  =  KHYPRG  +  1 
EM  =  EM  -  1. 

10  EN  =  EN  -  1. 

RETURN 

END 
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Figure  3-8 
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Random  number  generation  algorithm  for  hypergeometric 
distribution 
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3.8  GEOMETRIC  DISTRIBUTION 

k- 1 

The  geometric  distribution,  p.  =  p(l  -p)  ,  describes  ti.e 

number  of  trials  to  the  first  success  in  a  series  of  Bernoulli  trials.  For 
p  ^ .  25  ,  the  geometric  distribution  is  most  efficiently  sampled  by  a 
direct  solution  of  the  discrete  inverse  equation.  When  p<  .25,  it  becomes 
more  efficient  to  generate  a  geometric  variate  by  truncating  an  exponential 
random  number. 

Sample  Routines 

For  p  <  .  25: 


FUNCTION  KGEOM  (ALNQ) 
C  ALNQ  IS  -A LOG  (1-P) 

R  =  EXPRN  (R) 

KGEOM  =  1  +  R/ALNQ 

RETURN 

END 


For  p^  .25: 


FUNCTION  KGEOM  (P) 
A  =  P 
Q  =  1  -  P 
KGEOM  =  1 
R  =  UNFRN  (R) 

10  R  =  R  -  A 

IF  (R.  LT.  0)  RETURN 
KGEOM  =  KGEOM  +  1 
A  =A  *  Q 
GOTO  10 
END 
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3.9  PASCAL  OR  NEGATIVE  BINOMIAL  DISTRIBUTION 
The  Pascal  distribution, 

describes  the  number  of  successes  occurring  before  the  nth  failure  in  a 
series  of  Bernoulli  trials.  For  low  or  moderate  values  of  p  ,  the  Pascal 
distribution  is  efficiently  generated  by  a  direct  simulation  of  a  sequence 
of  Bernoulli  trials.  As  p  becomes  large  (p>.75),  it  becomes  more 
efficient  to  sample  by  generating  a  geometric  variate  for  the  number 
of  trials  to  each  of  the  n  failures. 

Sample  Routines 

For  p  <  .  75: 

FUNCTION  KPASCL  (P,  N) 

KPASCL  =  0 
DO  20  J  =  1,  N 
10  R  *  UNFRN  (R) 

IF  (R.  GT.  P)  GO  TO  20 
KPASCL  =  KPASCL  +  1 
GOTO  10 
20  CONTINUE 
RETURN 
END 

For  p> .  75: 

FUNCTION  KPASCL  (ALNP,  N) 

C  ALNP  IS  -ALOG(P) 

KPASCL  =  0 
DO  10  J  =  1,  N 
I  =  EXP  RN(R) /ALNP 
10  KPASCL  =  KPASCL  +  I 
RETURN 
END 
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Figure  3-11. 


Generate  R«-U(0, 1) 


Random  number  generation  algorithm 
for  Pascal  distribution 


Figure  3-12.  Random  number  generation  algorithm 
for  Pascal  distribution  (continued) 


3. 10  CAUCHY  DISTRIBUTION 


The  Cauchy  distribution, 


f(x)  =  - w-  ,  -  «  <  x  <« 

*[1  +  (x-m)2] 

represents  the  distribution  of  the  ratio  of  two  normally  distributed  numbers. 
It  also  represents  the  tangent  of  a  random  angle.  It  is  easily  generated  by  a 
rejection  technique  which  selects  x  and  y  uniformly  in  a  unit  circle,  then  cal¬ 
culates  the  tangent  x/y. 

Caution:  The  moments  of  the  Cauchy  distribution  are  infinite;  the  behavior 
of  Cauchy  variates  in  a  simulation  will  be  erratic. 


Sample  program: 

FUNCTION  COCHRN  (AMU) 

10  X  =  UNFRN(Y) 

Y  =  2.  ♦UNFRN  (X)  -  1. 

IF  (X  *  X  +  Y  *  Y.  GT.  1)  GO  TO  10 
COCHRN  =  AMU  +  Y/X 
RETURN 
END 


34 


3. 11  RAYLEIGH  DISTRIBUTION 


The  Rayleigh  distribution, 


f(x) 


» 


is  derived  as  the  radial  error  when  the  x  and  y  errors  are  independent  normal 
variates.  It  has  a  simple  inverse  which  provides  the  most  efficient  method  for 
generating  Rayleigh  variates. 

Sample  routine: 


FUNCTION  RAYLRN  (SIGMA) 

RAYLRN  =  SIGMA  *  SQRT  (2.  *EXPRN(R)) 

RETURN 

END 
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f(x)  - 


3. 12  GAMMA  DISTRIBUTION 


The  gamma  distribution 


f(x) 


-xx 

HtJ  e 
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describes  the  time  for  exactly  rj  events  to  occur  when  events  occur  at  a 
constant  rate  X  .  When  i?  is  an  integer,  there  is  a  simple  combination  tech¬ 
nique  for  generating  gamma  variates.  However,  as  the  gamma  distribution 
is  one  of  the  Pearson  family  of  distributions,  there  is  a  need  for  selecting 
gamma  variates  when  rj  is  non -integral  even  though  there  is  no  physical 
model  for  this.  This  is  a  much  harder  task  but  can  be  accomplished  by  a 
combination  of  the  usual  technique  for  the  integral  part  of  rj  with  a  composite 
rejection  technique  designed  to  select  from  x*e"x  where  f  is  the  fractional 
part  of  tj. 


Sample  routines: 

For  v  integer: 

FUNCTION  GAMRN  (ALAM,  NETA) 

Y  =  1 

DO  10  I  =  1,  NETA 
10  Y  =  Y  *  UNFRN  (Y) 

GAMRN  =  -  ALOG(Y)/ALAM 

RETURN 

END 

For  rj  general: 

FUNCTION  GAMRN(ALAM,  ETA) 

N  =  ETA 
F  =  ETA  -  N 
IF(F.EQ.O)  GO  TO  100 
10  R  =  UNFRN(R) 

IF  (R.LT.  F/(F  +  2.71828))  GO  TO  20 

Y  =  UNFRN(Y)  **  (1/F) 

IF  (UNFRN(R).  GT.  EXP(-Y))  GO  TO  10 
GO  TO  50 
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20 

50 
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80 
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Y  =  0 
GOTO  70 


Y  *  1.  +  EXPRNfY) 

IY(UNFRN(R).  GT.  Y**  (P-1.))  TO  TO  10 
IF(N.EQ.  0)  GO  TO  150 
Z  *1.0 


DO  80 1  «  1,  N 
z  *  Z*  UNFRN(Z) 
Y  =  Y  -  ALOGfZ) 
gamrn*  Y/ALAM 
return 


END 
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3. 13  BETA  DISTRIBUTION 


The  beta  distribution, 

,(x)  ’  e “h  (^rl) 

with  x  limited  to  the  interval  (a, b) ,  is  a  basic  statistical  distribution  fre¬ 
quently  encountered  for  bounded  variables.  The  parameters,  y  and  r\ , 
are  limited  to  positive  values.  Beta  variates  for  most  values  of  the  parame¬ 
ters  are  best  obtained  as  a  ratio  of  two  gamma  variiies.  If  y  and  rj  are 
both  small  integers,  a  beta  variate  may  also  be  generated  by  choosing 
y  +  rj  -  1  uniform  random  numbers,  arranging  them  in  order  of  increasing 
magnitude,  and  selecting  the  yth  random  number  as  the  beta  variate. 

Sample  routine 

FUNCTION  BETARN  (GAM,  ETA,  A,  B) 

Y  =  GAMRN  (1. ,  GAM) 

Z  =  GAMRN  (1. ,  ETA) 

BETARN  =  (Y/(Y  +  Z))  *  (B  -  A)  +  A 

RETURN 

END 
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Figure  3-16.  Random  number  generation  algorithm 
for  beta  distribution 


3. 14  PARETO  DISTRIBUTION 

The  Pareto  distribution,  f(x)  =  Xc*x’*”*,  has  a  simple  inverse 
which  provides  the  quickest  procedure  for  random  number  selection. 

Sample  routine 

FUNCTION  PRTORN  (EPS,  ALAM) 

PRTORN  =  EPS  ♦  UNFRN(R)  *  *  ( -1 .  /ALAM) 

RETURN 

END 
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3. 15  LOG-NORMAL  DISTRIBUTION 


The  log-normal  distribution 


f(x)  =  - -  exp  { -  [In  (x  -  c)  -  nf  j 

oj2*(x-c)  {  2 a  » 


describes  a  random  variable  whose  logarithm  is  normal.  It  is  a  simple 
matter  then  to  invert  this  transformation  to  generate  log-normal  variates. 

Sample  routine: 

FUNCTION  ALNMRN  (EPS,  AMU,  SIGMA) 

R  =  ANRMRN(R) 

ALNMRN  =  EPS  +  EXP  (SKHfA*R  +  AMU) 

RETURN 

END 
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Figure  3-18.  Random  number  generation  algorithm 
for  log-normal  distribution 


3.  16  FOLDED-NORMAL  DISTRIBUTION 
The  folded-normal  distribution, 

ft)  =  _»  [eW/*1  +  e  W/2o2]  , 

ojto! 

describee  the  distribution  of  the  absolute  value  of  a  normal  variate,  which 
provides  the  simplest  procedure  for  generating  from  the  distribution. 

Sample  routine 

FUNCTION  FNRMRN  (AMU,  SIGMA) 

FNRMRN  =  ABS  (AMU  4  SIGMA  *  ANRMRN(R)) 

RETURN 

END 
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Generate  R  *  N(0, 1) 


Figure  3-19.  Random  number  generation  algorithm 
for  folded -normal  distribution 


3.  17 


KODLDrS  DISTRIBUTION 


Kodlin  suggested  as  a  distribution  for  survival  time  data  the  functional 


form, 


f(x)  «  (fi  +  y  x)e'("*  +  1/2  y*i)  . 

This  Kodlin  form  has  a  moderately  simple  inverse,  and  thus  it  is  not  difficult 
to  generate  random  varities. 

Sample  routine 


FUNCTION  AKODRN  (ETA,  GAM) 

R  *  EXPRN  (R)  *  2.  *  GAM/ (ETA  **2) 
AKODRN  =  ETA/GAM  *  (SQRT(1.  +  R)  -  1.) 
RETURN 
END 


Figure  3-20.  Random  number  generation  algorithm 
for  Kodlin's  distribution 


3. 18  EXTREME  VALUE  DISTRIBUTIONS 


There  are  two  extreme  value  distributions.  The  first  is  for  the  maxi¬ 
mum  value, 


f«  [-  l(xy)  -  ]  , 

and  the  second  is  for  the  minimum  value, 

The  inverse  function  for  both  is  straightforward  and  provides  an  efficient 
selection  procedure. 

Sample1  routines 

For  the  maximum  value: 

FUNCTION  AMAXRN  (AMU,  SIG) 

R  =  EXPRN  (R) 

AMAXRN  =  AMU  -  SIG  *  ALOG  (R) 

RETURN 

END 

For  the  minimum  value: 

FUNCTION  AMINRN  (AMU,  SIG) 

R  -  EXPRN  (R) 

AMINRN  *  AMU  +  SIG  *  ALOG(R) 

RETURN 

END 
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Maximum  value:  f(x)  =  -  expL  i  (x-  n)  -  e”  ^x” 

o  [  ff 


START 


Generate  R «-  E(0, 1) 


=  -a  In  R  +  n 


Minimum  value:  f(x)  =  ^  exp|~  (x-  n)  -  e 


START 


Generate  R  «-  2(0, 1) 


=  or  In  R  +  m 


Figure  3-21.  Random  number  generation  algorithm 
for  extreme  value  distributions 


3. 19  WEIBULL  DISTRIBUTION 

The  Weibull  distribution,  f(x)  *  tj/x(x-e)^’1exp[-(x-^TVx  ],  iS  a  three' 
parameter  (c ,  X,  y)  family  of  empirical  distributions  having  wide  usefulness. 
Hie  random  variable  x  is  bounded  below  by  *  .  The  inverse  cumulative 
function  is  straightforward  and  provides  the  best  general  method  for 
generating  Weibull  random  numbers. 

Sample  routine: 

FUNCTION  WIBLRN  (EPS,  ALAM,  ETA) 

WDBLRN  =  EPS  +  (ALAM  *  EXPRN  (ALAM))  **  (1.  /ETA) 

RETURN 

END 
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Figure  3-22.  Random  number  generation  algorithm 
for  Weibull  distribution 


3. 20  JOHNSON  DISTRIBUTIONS 
3.20.1  Johnson  St.  Distribution 

,(x)= R+t*x'()]2}  * 

is  easily  generated  by  transforming  a  normal  variate.  (The  reverse  of  the 
transformation  used  in  deriving  this  Johnson  distribution.)  The  SL  dis¬ 
tribution  is  also  known  as  the  log-normal  (Section  3. 15). 

Sample  routine: 

FUNCTION  SLRN  (EPS,  GAM,  ETA) 

R  =  ANRMRN  (R) 

SLRN  =  EPS  +  EXP  ((R -GAM) /ETA) 

RETURN 

END 
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3.20.2  Johnson  Sp  Distribution 

The  Johnson  Sfi  distribution, 


,<x)  =  jki  ix-()  exp  }  ■ 

is  easily  generated  by  a  transformation  on  a  normal  variate. 
Sample  routine: 


FUNCTION  SBRN  (EPS,  ALAM,  GAM,  ETA) 
R  =  ANRMRN  (R) 

EX  =  EXP  ((R -GAM) /ETA) 

SBRN  =  EPS  +  ALAM  *  rX/(l.  +  EX) 

RETURN 

END 


I 


3.20.3  Johnson  Sg  Distribution 


Like  the  other  Johnson  family  distributions,  the  distribution. 


f(x)  = 


V(x+c)2  +  X2 


is  easily  selected  by  reversing  the  transform  which  generated  the  distribu 
tion  from  a  normal  distribution. 


Sample  program: 

FUNCTION  SURN  (EPS,  ALAM,  GAM,  ETA) 

R  =  ANRMRN(R) 

SURN  =  EPS  +  ALAM  *  SINH  ((R  -  GAM) /ETA) 

RETURN 

END 
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Figure  3-25.  Random  number  generation  algorithm 
for  Johnson  Sy  distribution 


3.21  PEARSON  DISTRIBUTIONS 


3.21.1  Pearson  Type  I  Distribution 

The  Type  I  distribution  of  the  Pearson  system  of  frequency  functions 
is  given  by 

f(a)  =  C(1  +  x/a.)ml  (1  -  x/a«»  )m2 

1  *  * 

where  C  is  a  normalization  constant.  The  limits  on  the  distribution  are 

-a,  <  x  <  a9  and  there  are  further  constraints  that  m1  >  -1  and  m9  >  *1. 

12  x  +  a-  1 

By  the  linear  transformation  Z  - - —  ,  the  Type  I  distribution  can  be 

*2+  l 

transformed  into  a  beta  distribution  which  may  be  derived  from  gamma  vari¬ 
ates  as  given  in  Section  3. 13. 

Sample  routine: 

FUNCTION  TYP1RN(EM1,  EM2,  Al,  A2) 

U  =  GAMRN  (1. ,  EM1+1. ) 

V  =GAMRN(1.,EM2+1.) 

TYP1RN  =  (Al  +  A2)*U/(U+V)  -  Al 

RETURN 

END 


i«  •  ca + x/»1)m  i  a  - 


Figure  3-26.  Random  number  generation  algorithm 
for  the  Pearson  Type  1  distribution 
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3.21.2  Pearson  Type  n  Distribution 

The  second  distribution  in  the  Pearson  family  is  given  by 

iw»cu-4>m  . 

a 


where  C  is  a  normalization  factor.  The  limits  on  the  distribution  are 
-a  <  x  <  a  and  m  >  -1.  This  is  a  special  case  of  Type  I  where  m^  = 
and  a^  =  As  such  it  may  also  be  derived  from  gamma  variates. 

Sample  routine: 

FUNCTION  TYPE2RN(EM,  A) 

U  =  GAMRN  (1. ,  EM+ 1) 

V  =  GAMRN  (1. ,  EM+1) 

TYP2RN  =  A*(U-V)/(U+V) 

RETURN 

END 


3.21.3  Pearson  Type  in  Distribution 

The  Pearson  Type  m  distribution  is  given  by  f(x)  =  C(1  +  x/a)  e“**, 
where  C  is  a  normalization  constant.  The  distribution  is  limited  to 
-a  <  x  <  a  (or  to  a  <  x  <  -a  if  a  is  negative)  and  is  further  constrained 
by  ya>  -1.  A  few  simple  transformations,  x  =  a(y-l)  and  X  =  a-y  ,  will 
turn  this  distribution  into  a  special  form  of  the  gamma  distribution 
f(y)  =C’yXe“Xy  . 

Sample  routine: 

FUNCTION  TYP3RN  (GAM,  A) 

P  =  GAM*A 
Y  -  GAMRN(P,P+1.) 

TYP3RN  =  A*(Y-1.) 

RETURN 

END 
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3.21.4  Pearson  Type  IV  Distribution 

The  Type  IV  distribution  of  the  Pearson  system  is  given  by 

f(*)=C(l  +  *VaVme-yt“'1(,t/a)  , 

where  C  is  a  normalization  constant.  By  a  trigonometric  transformation, 

* 

x  =  a  tan-1  (p -  tr/2),  the  function  can  be  transformed  into  f(p)  =  c’(sin o)Te~Y<0 
where  y  =  2m  -  2.  In  this  form  there  is  one  limit  on  the  parameters,  namely 
r  >  3,  while  p  ranges  from  0  to  v  .  Picking  from  this  function  can  be 
accomplished  by  a  selection  from  e"y<0,  truncated  at  p  =  it ,  followed  by  a 
rejection  conditioned  on  (sin  p)  . 

Sample  routine: 

FUNCTION  TYP4RN(EM, GAMMA,  A) 

DATA  PI/3. 1415962/HAFPI/1. 5707981/ 

R  =  2*EM-2 
10  PHI  =  EXPRN(R) 

PHI  =  AM0D(PHI/ GAMMA,  PI) 

IF  (UNFRN(R).  GT.  (SIN(PHI)**R))  GO  TO  10 
TYP4RN  =  A*TAN(PHI  -HAFPI) 

RETURN 

END 
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fW=C<l  +  x2/»V 


m  -y tan'Vx/a) 


1 


Figure  3-29.  Random  number  generation  algorithm 
for  the  Pearson  Type  IV  distribution 


3.21.5  Pearson  Type  V  Distribution 

The  fifth  typo  of  distribution  in  die  Pearson  system  of  frequency  func¬ 
tions  is  given  by  f(x)  =  C  x“pe"^x  ,  where  C  is  a  normalization  constant. 
The  range  of  the  argument  is  0  <  x  <  •  .  The  parameter  y  must  be  positive 
(for  y  <0,  -•  <x  <  0)  ,  and  p  must  be  greater  than  1 .  The  Type  V  random 
variable  x  is  the  inverse  of  a  gamma  variate;  this  provides  the  simplest 
means  of  picking  from  the  Type  V  distribution. 

Sample  routine: 

FUNCTION  TYP5RN  (P,  GAMMA) 

TYP5RN  =  l./GAMRN(GAMMA,P-l.) 

RETURN 

END 
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3.21.6  Pearson  Type  VI  Distribution 

Type  VI  of  the  Pearson  family  of  distributions  is  given  by 
f(x)  =  C(x-a)^  x  1  ,  where  C  is  a  normalization  factor  and  and  q2 
are  parameters  limited  by  q^  >  q2  +  1  >  0.  For  a  >  0  the  range  of  the 
distribution  is  a  <  x  <  •  while  for  negative  a  it  is  -®<  x  <  a.  By  the 
simple  transformation  x  =  a/y  the  distribution  is  converted  into  a  form  of 
the  beta  distribution 


(q«  -qo-2)  % 

f(y)=C’y  1  z  (1-y)  &  0  <  y  <  1 

which  can  be  obtained  from  two  gamma  variates  as  described  in  3. 13. 
Sample  routine: 

FUNCTION  TYP6RN(A,Q1,Q2) 

U  =  GAMRN(1.  ,Q1-Q2-1. ) 

V  =  GAMRN(1.  ,Q2+1.) 

TYP6RN  =  A*(U+V)/U 

RETURN 

END 
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3.21.7  Pearson  Type  vn  Distribution 

Type  VII  of  the  Pearson  family  of  distributions  is  given  by 

f(x)  =  C  (1  +x2/a2)"m, 

where  C  is  a  normalization  factor.  The  range  of  x  is  -  •  to  •  where 

m  must  be  greater  than  2.5.  By  setting  z  =  -j - g  the  distribution 

is  transformed  into  a  +  x 

gin)  =  C'  (1  -  z)'1/2 

which  is  a  special  case  of  a  Beta  distribution  with  y  =  m-1/2  and 

V  =  1/2.  The  beta  variate  z  can  be  obtained  as  a  ratio  of  two  gamma 

1/2  1/2 

variates,  z  =  u/(u+v)  .  As  x  =  a(l/z  -  1)  ,  we  have  x  =  a  (v/u) 

Now  v  is  a  gamma  variate  with  parameter  n  -  1/2.  This  special  case 
of  a  gamma  variate  can  be  obtained  from  v  =  y  /2  ,  where  y  is  a 
normalized  normal  variate.  This  gives  x  =  ay  (l/2u)  7  .  Selection 

from  the  Pearson  Type  VII  is  achieved  by  combining  the  above  transform¬ 
ations  with  the  selection  routines  for  the  gamma  and  normal  variates. 

Sample  Routine 

FUNCTION  TYP7RN(A,  EM) 

Y  =  ANRMRN(Y) 

U  =  GAMRN  (.5,  EM -.5) 

TYP7RN  =  A*Y/SQRT(U) 

RETURN 

END 
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3.21.8  Type  Vm  Pearson  Distribution 

The  eighth  distribution  in  the  Pearson  family  is  given  by 

f  M  =  C  (1  +  x/a )•” 

where  C  is  a  normalisation  constant.  The  range  of  x  is  -a  <  x  <  0 
(or  0  <  x  <  -a  for  a  negative)  while  the  range  of  m  is  Os  ms  1. 
ft  we  set  y  =  (1  +  x/a) ,  the  distribution  becomes 

i(y)  =  C'  y”m  where  0<  y<  1. 


This  form  of  the  distribution  has  a  simple  inverse. 
Sample  Routine 

FUNCTION  TYP8RNCA,  EM) 

R  =  UNFRN(R) 

TYP8RN  =  A*(R**(1./(1.EM)) -1.) 

RETURN 

END 
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f  (*)  =  C  (1  +  x/a)“m 


Random  number  generation  algr  rithm 
for  the  Pearson  Type  VTn  distribution 


3.21.9  Pearson  Type  IX  Distribution 

The  Pearson  Type  DC  distribution  is  given  by 

m  =  C  (1  +  x/a)m, 


where  C  is  the  normalization  factor.  The  range  of  x  is  -a  to  0  while 
m  must  be  greater  than  zero.  This  function  has  a  simple  inverse. 

Sample  Routine 


FUNCTION  TYPE9RN(A,  EM) 

R  =  UNFRNCR) 

TYP9RN  =  A*(R**  (1.  /(EM  ♦  l.))-l.) 

RETURN 

END 


f- 
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i  (*)  =  C  (1  +x/a)m 


Figure  3-34.  Random  number  generation  algorithm 
for  the  Pearson  Type  IX  distribution 
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3.21.10  Pearson  Type  X  Distribution 

The  Pearson  Type  X  distribution  is  a  form  of  die  exponential 
distribution  given  by 

f(a)  =  1/e  e”x ^  ;  x  *  0 

This  is  easily  obtained  from  the  standard  exponential  distribution 
routines. 


FUNCTION  TP10RNC3IGMA) 
TP10RN  =  SIGMA*  EXPRN  C3IGMA) 
RETURN 
END 


Sample  Routine 


f (x)  =  l/o  e  "x/a 


Figure  3-35.  Random  number  generation  algorithm 
for  the  Pearson  Type  X  distribution 
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3.21.11  Pearson  Type  XI  Distribution 

The  eleventh  in  the  series  of  Pearson  distribution  is  given  by 

f(x)  =  C(b/x)m 


where  C  is  a  normalization  factor.  The  range  of  x  is  limited  to 
b  <x  <•.  The  parameter  m  is  greater  than  1  .  This  distribution  has 
a  simple  inverse. 

Sample  Routine 

FUNCTION  TP11RN(B,  EM) 

R  =UNFRN(R) 

Y  =R**(1./(EM-1.)) 

TP11RN  =  B/Y 

RETURN 

END 
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f(x)  .  C(b/x)m 


Figure  3-36.  Random  number  generation  algorithm 
for  the  Pearson  Type  XI  distribution 


3.21.12  Pearson  Type  XII  Distribution 

Type  XII  of  the  Pearson  system  of  distributions  is  given  by 


J°(ysx+/r)+*i('/v<^) 

where  C  is  a  normalization  factor,  o  is  the  standard  deviation,  and 
2  3 

fj  =  M3/M2  (skewness).  The  range  of  x  is 

*  o ) <  x  <  a{J 3+/J j  ■  . 

By  setting 

m  =  , 

a  s  +  ^3j)  f  and 

b  =  -  J 8 j) » 

the  distribution  becomes 


«  ■  c  fef  • 

By  setting 


the  distribution  transforms  to  f(y)  =  C’y1*1  (l-y)’m  which  is  a  special  case 
of  the  Beta  distribution. 
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Stmplrgcgttg* 


FUNCTION  TY12RN(SIGMA,  BETA1) 
R  *  SQRT(BETAl) 

S  -  SQRT(BETAl+3) 

EM  *  R/S 
A  =  SKaiA*(R+S) 

B  =  SIGMA*  (8 -R) 

Y  *  BETARN(EM+1, 1-EM) 

TP12RN  *  (B+A)*Y-A 

RETURN 

END 


3.22  HISTOGRAM  DISTRIBUTIONS 


Frequently,  empirical  data  regarding  a  probability  distribution 
is  obtained  in  a  histogram  form.  That  is,  intervals  (xQ,  xj)  ,  (x^,  x^), .... 

(x  .,  x  )  and  probabilities  p.,p0 . ,  p,  are  given  such  that  p.  is 

the  probability  that  the  variable  x  is  found  in  the  interval  from  Xj  ^  to 

x. .  (It  is  presumed  that  the  histogram  is  normalized,  i.e.  1  p  =  1.) 

1  i=l  1 

Within  each  interval  it  is  assumed  that  the  probability  is  constant. 

Selecting  a  random  number  from  such  a  histogram  distribution  is 
simple.  It  is  necessary  first  to  select  the  interval  in  which  the  random 
number  falls,  and  then  to  choose  where  in  that  interval  the  random  nv.mber 
lies.  This  is  basically  an  inverse  distribution  technique.  Selection  of 
the  interval  i  is  accomplished  by  generating  a  uniform  random  number  and 
subtracting  off  successive  values  of  p^  .  The  value  of  i  when  this  result 
first  goes  negative  is  the  desired  interval  index.  Generation  of  a  second 
uniform  random  number  and  scaling  it  to  fit  in  the  interval  from  Xj_^  to  x. 
completes  the  task. 

A  more  efficient  (much  more  efficient  if  the  size  of  the  data  table 
is  large)  generator  can  be  produced  if  it  is  possible  to  cast  the  histogram 
<ata  in  a  form  such  that  P^  =  Pg  =  •  •  •  =  Pn  =  1/n  by  choosing  values  of 
x.  appropriately.  Such  a  representation  is  known  as  equal  probability 
bins.  This  greatly  simplifies  selection  of  the  interval  i  as  all  n  intervals 
have  the  same  probability.  Successive  subtraction  of  values  of  p.  is 
no  longer  needed  and  can  be  replaced  by  a  direct  calculation  of  i  from 
a  uniform  random  number. 

In  the  sample  Fortran  routines  below,  the  array  Xfl)  is  presumed 
to  contain:  X(l)  =  xQ  ,  X(2)  =  x^  ,  . . . . ,  X(N  +  1)  =  xn»  In  the  first  routine 
use  is  made  of  the  fact  that,  at  the  conclusion  of  selection  of  i,  R  will  be 
uniformly  distributed  between  0  and  -  pi  . 
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Sample  Routines 

For  general  histogram  selection 

FUNCTION  HISTRN  (N,  X,  P) 

DIMENSION  X(N),  P(N) 

R  =  UNFRN  (R) 

DO  10  I  =  1,  N 
R  =  R  -  P(I) 

IF  (P  .  LT  .  6)  GO  to  20 
10  CONTINUE 

20  HSTRN  =  X(I)  -  R  *  (X  (i  +  1)  -  X  (I))/P(I) 
RETURN 
END 

For  selection  with  an  equal  probability  bin  histogram 

FUNCTION  HSTRN  (N,  X) 

DIMENSION  X(N) 

R  =  N  *  UNFRN  (R)  +  1 
I  =  R 
R  =  R  -  I 

HSTRN  =  X(t)  +  R  *  (X(I  +  1)  -  X(l)) 

RETURN 

END 
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GENERAL  TECHNIQUES  FOR  GENERATING  RANDOM 
NUMBERS  FROM  DESIRED  DISTRIBUTIONS 

When  given  a  particular  distribution,  f(x),  and  the  task  of 
selecting  random  numbers  distributed  according  to  that  function,  the 
investigator  has  a  large  number  of  possible  alternatives  at  his  disposal. 

The  primary  task  is  to  derive  a  method  which  will  accomplish  the 
desired  selection.  A  secondary  task  is  to  choose  the  method  which  is 
least  time-consuming  computationally. 

Unfortunately,  it  is  not  possible  to  give  a  straightforward 
methodology  for  deriving  random  number  generation  techniques  which 
can  be  applied  in  all  or  even  in  most  cases.  The  situation  closely 
parallels  that  of  finding  an  integral  of  an  arbitrary  function.  When  one 
encounters  the  need  to  integrate  an  unfamiliar  function,  the  first  step,  of 
course,  is  to  try  to  look  it  up  in  a  table  of  integrals.  That  failing, one  must 
try  to  simplify,  transform  variables,  integrate  by  parts,  use  trigonometric 
substitutions,  or  employ  other  similar  tricks  to  reduce  the  integral  to  a  familiar 
form.  There  is  no  guarantee  of  success,  and  much  depends  on  the  ingenuity 
and  experience  of  the  researcher.  When  all  else  fails  you  can  "grind  out" 
a  numerical  solution. 

Faced  with  the  task  of  generating  random  numbers  from  an  unfamiliar 
distribution,  a  similar  procedure  is  needed.  The  first  step  is  to  try  to  look 
it  up  somewhere  —  such  as  in  Section  3  of  this  report.  If  not  found  there, 
there  are  a  number  of  techniques  —  inverse,  rejection,  transformations, 
combinations,  etc.  available.  These  are  described  in  this  Appendix.  There 
is  no  guarantee  of  success  in  using  them,  and  the  experience  and  ingenuity 
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of  the  analyst  is  very  important.  As  a  final  resort,  there  are  numerical 
methods  which  can  be  applied. 

The  following  description  of  general  techniques,  while  not  universally 
applicable  should  give  the  reader  some  notion  of  how  to  proceed  in  deriving 
random  number  generation  algorithms. 

A.  1  THE  INVERSE  METHOD 

The  first  technique  which  one  should  consider  is  the  inverse.  To 
apply  the  inverse  method,  the  distribution  function  is  integrated  to  give 
the  cumulative  distribution,  F(x)  =  Jx  f(x’)dx\  This  is  the  probability  of 
selecting  a  number  less  than  or  equal  to  x.  This  is  equated  to  the  proba¬ 
bility  of  selecting  a  random  number,  R,  from  the  uniform  distribution. 
Thus,  F(x)  =  fX  f(x')  dx'  =  R.  The  question  then  is  whether  or  not  this 

-*  _i 

equation  has  a  simple  closed -form  solution,  x  =  F  (R).  If  the  inverse 
function  exists,  then  it  is  a  solution  to  our  task,  for,  if  R  is  distributed 
uniformly,  then  x  =  F-1(R)  is  distributed  according  to  f(x).  If  F_1(R) 
not  only  exists, but  is  also  moderately  simple  to  compute,  it  is  most  likely 
the  most  efficient  way  to  generate  the  desired  random  numbers. 

A.  2  REJECTION  TECHNIQUE^ 

If  the  inverse  function  cannot  be  easily  calculated,  then  the  rejection 
technique  should  be  considered.  Suppose  that  the  function,  f(x),  has  a 
maximum  value  M  where  x  varies  over  the  range  of  interest  from  a  to  b. 
Random  numbers  are  then  chosen  by  the  following  two-step  procedure. 

•  Select  x  from  a  uniform  distribution  on  the  interval  (a,b) 

•  Select  a  second  uniform  random  number,  y,  and  accept 
the  value  x  only  if  y  <[f(x)]/M. 

If  x  is  rejected,  then  go  back  to  the  first  step  to  select  a  new  x  and  con¬ 
tinue  this  procedure  until  some  value  of  x  is  accepted.  The  probability  of 
selecting  x  in  the  first  step  is  f l/(b-a)  ]  dx,  while  the  probability  of 
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acceptance  at  the  second  step  is  f(x)/M.  Thus  the  x  values  will  be  genera¬ 
ted  with  the  desired  probability  f(x)  dx. 

The  constant  term  l/[M(b-a)]  represents  the  efficiency  of  the  rejec¬ 
tion.  Its  reciprocal,  M(b-a),  is  the  average  number  of  trials  the  rejection 
technique  will  require  to  generate  a  single  random  number  and  is>  therefore* 
linearly  proportional  to  the  computation  time  required.  If  M(b-a)  is  very 
large,  the  rejection  technique  is  too  inefficient  and  a  better  technique  should 
be  sought. 

The  rejection  technique  need  not  be  based  on  variables  from  a  uniform 
distribution  but  can  be  developed  from  other  distributions.  For  example 
the  fact  that 

-x2/2  1/2  -x 

e  <  e  •  e 


can  be  used  to  develop  a  rejection  technique  for  picking  from  a  normal  distribution. 
First  select  x  from  the  exponential  distribution  e"x.  Then  accept  x  if 
a  second  (uniform)  random  number 


e'*2/2 

y  <~ - II 

Je  •  e  x 


e-(*-i)2/2 


The  essential  ingredient  of  the  rejection  technique  is  to  find  a  second  dis¬ 
tribution  function,  g(x),  for  which  a  selection  procedure  is  known  and  such 
that  f(x)  <  C  g(x).  Selection  of  x  from  g(x)  is  followed  by  acceptance  if 


y  <rira  • 

The  average  number  of  trials  needed  for  an  acceptance  is  C.  Note  that  if 
g(x)  is  close  to  f(x),  then  C  will  be  close  1  and  the  technique  will  be  very 
efficient. 


A.  3  TRANSFORMATION 

To  simplify  the  derivation  of  inverse  or  rejection  methods,  it  is  best 
to  transform  the  random  variable  into  its  simplest  form.  Thus,  if  one  had 
f(x)  =  g(Xx  +  c),  one  would  first  make  the  substitution,  y  =  Xx  +  e,  then 
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search  for  a  technique  for  generating  numbers  from  g(y).  After  generating 
a  random  number  y,  set  x  -  (y- c  )/A  to  get  the  desired  random  variable. 

In  doing  transformations  correctly  we  must  be  careful  to  transform  not  just 
the  function  f(x)  but  the  probability  f(x)  dx.  Thus,  properly,  we  have 
f(x)  dx  =  g(Xx  +  c)  dx  =  g(y)  dx  =  g(y)  dyA  as  the  substitution  y  =  Xx  +  c 
implies  dy  =  X  dx.  The  correct  normalized  distribution  for  y  is  then 
1/X  g(y).  As  a  second  example,  assume  f(x)  dx  =  2x  dx.  Try 
the  transformation  y  =  x  .  As  dy  -  ?x  dx,  f(x)  dx  =  2x  e"x  dx  =  e"y  dy. 
Therefore,  selecting  y  from  the  exponential  e  ^  and  taking  x  =  Jy  will 
give  a  random  x  from  f(x). 

A.  4  COMBINATION  OF  RANDOM  VARIABLES  ^ 

As  a  step  beyond  transformations,  consider  various  combinations 
of  random  variables  such  as  adding  subtracting,  or  multiplying  two 
random  numbers,  taking  the  maximum  or  minimum  of  several  random 
numbers,  etc.  The  results  of  such  combinations  follow  no  intuitive  pattern 
but  must  be  worked  out  through  the  laws  of  probability.  For  example,  the 
sum  of  two  uniform  random  numbers  has  a  triangular  distribution, 
f(x)  =  1  -  |x  -  1|  while  the  product  has  the  distribution,  f(x)  =  -  In  x. 

More  complex  examples  seem  even  farther  removed  from  simple  ration¬ 
ality.  If  x  and  y  are  random  numbers  from  the  gamma  distributions, 
l/r(n)  x11”1  e~x  and  l/r(m)  ym_1  e"x,  then  z  =  x/(x+  y)  has  a  beta 
distribution  T(m  +  n)/r(m)r(n)  zn“*  (1  -  z)m"*.  However,  the  beta 
distribution  may  also  be  obtained  by  taking  n  +  m  -  1  uniform  random 
numbers,  arranging  them  in  increasing  order,  and  selecting  the  nth  num¬ 
ber  in  the  sequence.  Thus,  although  combinations  can  be  a  very  powerful 
method  for  transforming  simple  random  variables  into  selections  from 
other  distributions,  it  is  impossible  to  give  guidelines  or  to  arrive  at  a 
methodology  for  determining  the  proper  combination  needed  to  arrive  at  a 
desired  distribution.  The  investigator  must  simply  learn  the  frequently 
used  combinations  and  must  use  his  inventiveness  when  confronted  with  an 
unfamiliar  distribution. 
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A. 5  COMPOSITION  TECHNIQUE^®* 

Another  method  of  general  applicability  is  the  composition  technique.  If 
the  desired  distribution  can  be  written  as  a  (generalized)  integral  over  a 
family  of  density  functions,  then  the  sampling  can  be  accomplished  in  a  two- 
stage  process.  On  the  first  step,  a  particular  density  function  is  selected 
from  the  family,  and  on  the  second  step,  the  desired  random  number  is 
drawn  from  the  particular  density  function.  In  the  usual  application  of 
this  technique,  the  desired  distribution  is  broken  down  into  discrete  parts, 
generally  on  separate  intervals. 

A.  6  NUMERICAL  METHODS 

If  no  exact  method  can  be  derived,  there  is  a  numerical  technique 
which  can  be  used.  This  consists  of  generating  the  cumulative  function, 
solving  for  its  inverse  numerically,  tabulating  the  inverse,  and  then  gener¬ 
ating  the  random  numbers  from  the  tabulated  data.  If  equal  probability  intervals 
are  used  in  tabulating  the  inverse,  then  generation  from  the  tabulated  data 
can  be  quite  fast.  It  does,  however,  require  a  certain  amount  of  computer 
storage  to  hold  the  tabulation. 

Improvements  in  the  accuracy  of  numerical  inverses  can  be  made  by 
using  Chebyshev  interpolating  polynomials/6*  For  some  functions  with  long 
tails,  the  tabulated  inverse  must  be  replaced  with  some  sort  of  approximating 
function  in  the  tail  of  the  distribution  to  achieve  reasonable  accuracy. 

A.  7  MARSAGLIA  TECHNIQUE*3-5* 

If  a  particular  distribution  is  very  central  to  a  frequently  used  simu¬ 
lation  program  and  the  generation  subroutine  will  be  called  a  great  many 
times  to  produce  random  numbers,  it  may  be  worthwhile  to  design  a  very 
fast  selection  procedure  to  reduce  the  computer  time  needed.  A  number  of 
super -efficient  techniques  have  been  developed  by  G.  Marsaglia.*3*  These 
are  based  on  composition  methods  where  the  function  is  expressed  as 
the  sum  of  three  or  more  parts.  The  parts  having  highest  probability  are 
fast  to  select  from  and  the  parts  difficult  or  slow  to  select  from  have  very 
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small  probability.  In  one  of  Marsaglia's  methods,  the  function  is  broken 
into: 

•  A  histogram 

•  A  collection  of  saw-toothed  functions  where  an  efficient 
rejection  technique  selects  from  the  'almost -linear’  dis¬ 
tribution  of  each  sawtooth. 

•  The  tail  of  the  distribution. 

This  method  is  very  fast  but  requires  moderate  amounts  of  computer  storage. 
In  another  method  distributions  are  fitted  to  an  approximation  of  the  form 
C(M  +  Uj  +  Ug  +  Ug),  where  M  is  a  discrete  variable  and  the  u's  are  uniform 
variables.  A  small  fraction  of  the  time  a  more  lengthy  rejection  procedure 
is  needed  to  correct  the  error  in  the  approximation.  This  method  is  fairly 
fast  without  great  storage  requirements. 

These  methods  have  been  applied  very  successfully  to  the  exponential 
and  normal  distributions.  They  do,  however,  require  considerable  effort 
in  manhours  to  develop  and  thus  should  be  applied  to  other  distributions  only 
when  the  payoff  can  justify  it. 
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MIRAN  -  A  MACHINE  INDEPENDENT  PACKAGE  FOR  GENERATING 
'  UNIFORM  RANDOM  NUMBERS 

B.  1  GENERAL  DISCUSSION 

The  standard  technique  for  producing  uniform  random  numbers  on 
modern  high-speed  computers  is  an  algorithm  known  as  the  multiplicative 
congruential  method.  This  method  is  expressed  mathematically  as 

Rn+1  =  X-Rr  (modulo  P)  . 

Since  the  R’s  are  integers  ranging  from  1  to  P-1,  successive  real  random 
numbers  uniformly  distributed  from  0  to  l  are  generated  by  dividing  Rn  by  P. 
The  properties  of  this  technique  as  a  random  number  generator  (RNG)  are 
highly  dependent  on  the  choice  of  the  generator,  X,  and  the  modulus,  P. 
Unfortunately,  there  are  many  RNGs  in  current  use  which  do  not  approximate 
randomness  closely  enough  to  be  sufficient  for  all  Monte  Carlo  calculations 
and,  what  is  far  worse,  do  manage  to  pacs  some  of  the  simple  tests  for 
randomness.  There  are,  however,  several  choices  of  X  and  P  which  have 
been  thoroughly  tested,  both  theoretically^  and  through  many  years  of  actual 
use  in  Monte  Carlo  calculations,  and  which  appear  to  be  sufficiently  random 
for  general  usage. 

For  reasons  of  convenience  and  efficiency,  P  is  generally  taken  to 

be  2m  where  m  is  the  number  of  bits,  excluding  the  sign  bit,  in  a  single 

word  on  the  particular  computer  being  used.  The  generation  process  starts 

with  a  fixed  generator,  X,  and  a  starting  value,  Rq.  The  full  product 

from  the  multiplication  of  X  and  Rq  would  usually  fill  two  computer  words; 

however,  the  modulo  P  in  the  algorithm  means  that  we  only  need  the  single 

word,  R, ,  comprising  the  low  order  half  of  the  X  •  R  product.  The  random 
l  o 

number  gen3ration  is  completed  by  converting  R^  to  a  real  variable  and 
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dividing  by  P.  replaces  RQ  in  storage  in  the  random  number  subroutine 
and  the  process  is  ready  to  begin  anew. 

In  this  sort  of  a  process  there  have  been  two  barriers  to  developing 
a  Fortran  RNG  subroutine  which  would  be  independent  of  the  particular  com¬ 
puter  for  which  it  was  designed.  The  first  is  the  modulus  P,  which  varies 
from  computer  to  computer  as  the  word  length  varies.  [Choosing  a  universal 
value  of  P  to  fit  the  smallest  computer  is  not  a  good  solution  as  the  proper¬ 
ties  of  a  RNG  become  less  random  as  P  is  made  smaller,  to  the  extent  that 
Coveyou  and  MacPherson^  consider  them  questionable  for  P  =  2^ 

*JR 

(IBM  360  series)  and  borderline  for  P  =  2  (IBM  7090,  Univac  1108,  etc.).  ] 
The  second  problem  is  that  the  sign  bit  of  may  need  to  be  cleared  follow¬ 
ing  the  multiplication.  Clearing  the  sign  bit  generally  requires  some  trickery 
in  Fortran  which  varies  from  computer  to  computer  as  the  mode  of  represen¬ 
tation  (one's  complement,  two's  complement,  uncomplemented,  etc.)  of 
negative  n  imbers  varies. 

The  way  around  these  obstacles  is  to  use  an  explicit  multiple  pre¬ 
cision  representation.  The  integers  and  operations  involved  in  the  RNG 
algorithm  are  separated  into  component  parts  in  such  a  way  that  all  operations 
are  kept  within  a  single  computer  word  and  no  overflows  into  the  sign  bit  are 
made,  thus  avoiding  the  sign-clearing  problem.  Through  multiple  precision 
a  sufficiently  large  modulus  for  good  RNG  properties  may  be  used  even 
though  the  actual  computer  word  size  is  small.  An  initialization  call  must 
be  made  to  convey  to  the  RNG  the  maximum  integer  allowed  on  the  particular 
computer  being  used  so  that  it  can  set  up  an  appropriate  multiple  precision 
representation. 

The  advantage  of  a  RNG  that  is  machine  independent  is  simple:  it 
greatly  facilitates  the  exchange  and  checkout  of  Monte  Carlo  programs  between 
different  computers.  The  price  paid  for  this  advantage  is  also  simple:  it 
is  a  much  slower  method  of  producing  random  numbers.  However,  it  is 
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still  fast  enough  (several  thousand  random  numbers  generated  in  one  second) 
that  the  time  difference  will  not  be  noticed  in  most  Monte  Carlo  applications. 

B.  2  CHOICE  OF  A  SPECIFIC  ALGORITHM  FOR  MIRAN 

The  work  of  Coveyou  and  MacPherson^  has  provided  a  thorough 
theoretical  analysis  of  many  commonly  used  RNGs.  They  show  that  the  cor¬ 
relation  properties  of  a  RNG  are  strongly  dependent  on  the  modulus  P. 

31  35 

For  values  of  P  =  2  or  2  ,  there  must  necessarily  be  a  waviness  or 
graininess  to  the  joint  distribution  of  two,  three,  and  four  consecutive  ran¬ 
dom  numbers  that  could  lead  to  incorrect  results  for  some  Monte  Carlo  cal- 

47 

culations  For  P  =  2  ,  the  departures  from  true  randomness  are  small 

enough  as  to  be  negligible  for  practical  calculations.  Among  the  specific 

15 

generators,  X  ,  tested  by  Coveyou  and  MacPherson,  there  is  one,  X  =  5  , 

which  has  good  statistical  properties  and  which  may  be  easily  produced  by 
a  machine  independent  subroutine.  (In  a  subroutine  designed  for  use  on  com¬ 
puters  of  varying  word  length,  specifying  a  fixed  47 -bit  integer  through 

15 

data  statements  would  be  difficult.  However,  5  may  easily  be  produced 

by  multiplying  5  ’s  after  the  exact  multiple  precision  representation  needed 

47  15 

has  been  established.)  In  addition  the  choice  of  P  =  2  and  X  =  5  has 
an  added  advantage:  this  particular  choice  of  a  RNG  has  seen  long  usage 
(several  thousand  hours  on  a  CDC  1604  at  Oak  Ridge  National  Laboratory) 
in  Monte  Carlo  computations  without  any  apparent  problems. 

B.  3  MULTIPLE  PRECISION  REPRESENTATION 

In  the  basic  algorithm  used  by  MIRAN.  X  and  the  R  values  will 

n 

be  47 -bit  integers.  This  may  exceed  machine  capacity.  To  keep  all  arith¬ 
metic  operations  from  overflowing  a  single  machine  word,  these  integers 
are  stored  in  an  array  wherein  each  word  of  the  array  constitutes  a  ’digit' 
in  a  representation  of  the  integer  to  a  particular  base.  This  basis,  called 
BASE,  is  chosen  at  execution  time  so  that  (BASE)  does  not  exceed  the  maxi¬ 
mum  integer  allowed  on  the  particular  computer  being  used.  Thus,  for 
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example,  on  a  machine  with  35 -bit  words  (unsigned),  BASE  would  be  2 
and  each  47-bit  integer  would  be  broken  down  into  3  words  as  follows: 


17 


47 -bit  Integer  Multiple  Precision  Representation 


blb2*  *  * 

*  *b13b14 - b30b31 - b47 

+  0, . . 

*  *0  bi - bi3 

word  3 

+  0. . . 

'  *0  bi4 - b30 

word  2 

+  0.  •  • 

•  •  0  bjj#  • « • 

word  1 

Note  that  the  'digits'  are  stored  in  the  array  in  'reverse'  order,  i.  e. , 
word  1  is  the  least  significant  17  bits  of  the  number.  Also,  since  17  does 
not  go  evenly  into  47,  the  last  word  contains  only  13  bits. 

Arithmetic  in  a  multiple  precision  representation  is  carried  out  in 
the  same  manner  as  arithmetic  is  normally  done  by  hand.  The  addition  of 
two  numbers,  for  example,  is  done  digit  by  digit.  When  two  'digits',  or  words, 
are  added  there  may  be  an  overflow  into  the  18th  bit  of  the  result.  This  must 
be  detected,  the  overflow  cleared  out,  and  a  carry  of  1  added  into  the  next 
higher  'digit'.  Multiplication  is  slightly  more  complex.  It  is  again  carried 
out  digit  by  digit  and  the  resulting  products  are  added,  keeping  them  in  appro¬ 
priate  columns,  to  get  the  final  product.  The  multiplication  of  two  'digits’ 
produces,  of  course,  a  two-digit  product  which  is  initially  contained  in  a 
single  computer  word.  This  must  be  broken  down  into  a  high-order  digit  and 
a  low-order  digit  with  the  high -order  digit  being  added  into  the  next  higher 
column  of  the  result.  As  each  column  is  added,  a  carry  over  into  the  next 
higher  column  may  be  needed.  Thus,  in  our  example  where  three  words  were 
used  for  each  integer,  nine  multiplies  and  several  additions  would  be  needed 
to  form  the  six -word  full  product  as  schematized  below. 
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'33 


d3 

d2 

dl 

d3 

di 

di 

hll 

Si 

h21 

l21 

h31 

*31 

h12 

*12 

h22 

*22 

h32 

*32 

h13 

*13 

h23 

*23 

*33 

S5 

S4 

S3 

S2 

S1 

6 


where  and  are  the  high  and  low  order  parts  of  *he  product  of 
dj  and  dj . 

B.  4  USE  OF  MIRAN  PACKAGE 
Initialization: 

Before  generating  any  random  numbers,  it  is  necessary  to  make  an 
initialization  call.  This  is  done  by  the  statement 

CALL  RANSET  (MAXINT,  NSTART) 

where  MAXINT  is  the  maximum  integer  allowed  on  the  computer  (or  compiler) 
being  used.  NSTART  is  the  starting  value,  Rq,  to  be  used  in  the  random 
number  sequence.  If  NSTART  is  less  than  or  equal  to  0 ,  a  default  value 
of  2001  is  supplied  for  NSTART.  If  NSTART  is  even,  the  next  higher  odd 
number  will  be  used. 
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For  example  MAXINT  =  2-1  on  a  1108,  2  0  - 1  on  a  CDC-6600,  etc. 
Good  values  for  NSTART  are  any  cdd  integer  although  frequent  use  of 
small  odd  integers  is  not  recommended  for  calculations  employing  a  re¬ 
latively  small  number  of  random  numbers. 

The  random  numbers  are  generated  in  subroutine  URAND  which  may 
be  used  as  either  a  function  subroutine  or  as  an  ordinary  subroutine  return¬ 
ing  a  value.  Thus,  either 

CALL  URAND(R) 
or 

R  =  URAND  (X) 

will  store  a  uniform  random  number  in  R.  (Note  that  in  the  second  form 
the  same  random  number  will  also  be  stored  in  X.  Thus,  X  must  be  a 
Fortran  variable  and  not  a  constant.) 

Limitations  of  MIRAN: 

MIRAN  will  work  on  all  computers  where  MAXINT  is  greater  than 
94 

1023  and  less  than  2  .  (These  limits  are  practical  and  not  theoretical  and 

could  be  extended  if  it  were  ever  necessary. ) 

3.  5  MIRAN  PROGRAM  DETAILS 

The  Fortran  listings  of  the  two  MIRAN  routines  URAND  and  RANSET 
are  presented  in  Figures  B-l  and  B-2.  The  accompanying  logic  flow  is  de¬ 
tailed  in  Figures  B-3  and  B-4.  Additional  explanation  of  the  last  step  in  the 
URAND  logic  is  provided  below. 

The  two  subroutines  URAND  and  RANSET  communicate  through  a 
labelled  common,  MIRNG  which  contains 

RAN(10)  -  An  array  containing  the  'digits'  of  the  current  (or  last) 

multiple  precision  random  integer 
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RtAL  FUNCTION  UnANDCFRAN) 

COMMON  /MJRNG/  KAN(10)#GFN{ 10)#NWRD#BaSE*MO0*F*ASE,FMUD 
DIMENSION  SUM(IO) 

INTEGER  RAN, G£N, RASE » CARRY* SUM* PROD* HPRQO 
00  30  I8al *N#HD 
RuM(I8)«0 
00  1  iGal,MNRu 

N2«NWR0-lG*i 
00  1  IRal,N2 
TSaIRMG-1 

»HCJ«nANUR)*GEn(IG) 

HPRU0«PH00/BA5E 
LPR0D»PHOd-mPN0U*BA8F 
S JM ( IS) *3UM ( I S ) ♦•  HROO 

TF  (IS.Ll.N^Ro)  SUM(I8*l)«PuM(I8tl)THPROO 
1  CONTINUE 
N2«NWR0-1 
00  5  IR»1,N2 
CARRYi80M(I8)/8ASE 
8UM(I8)«SjMII8)-CARHV*RASE 
8uMci8ti)«SuM(i8*i) ♦Carry 
5  CONI INUt 

SUM(NNK0)a80H(NMR0)«M0O*(8UM(NMKD)/M0D) 

00  iO  I8»l *NwRO 
20  PANCI8)a8uM(I3) 

FnANa8UM(i) 

00  10  1 3*2  *  N*RO 

10  FNANbF-»AN/F0AoR  +  9UM(13J 

pRANbFRAN/FmOd 
l)N  ANObFRAn 

return 

End 


Figure  B-l.  Fortran  listing  of  URAND 
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8URR0UTJNE  RAN3ET(MAXINT»N8TRT) 

COMMON  /MXRnG/  RAN(10),UFN(10),NWRD,BASEfMOD.FBA8E»FMOO 
INTEGER  RAN, G*.N, BASE, CARRY, KEH 

M4XI«MAXInT/U 

IH»0 

BA3t»l 

9 9  IF  i«A3fc.i.T,MAXl)  Gu  TO  10ft 
BASe»8ASF*A 

Gu  m  99 

100  BA3Ea*«*I* 

FBAbEaPAlE 

NwRQa  4  //IB* 1 

OEMaa7-IB*(M*iKO-l) 

MUOa?**RFM 

fmho.moo 

Db  101  N* 1 , 1  0 
RANfNJaO 

101  GLN(N)aO 
GtNll Ja5 

Du  *oo  i ■ s  *  i a 
CARRYsO 

Du  190  N« 1 , N*WO 

GE^(N)aCFNCM*S4CAR«Y 

CARRYaO 

IF  CDEM(N).LT.BA8E)  GO  TO  190 
CARRVaGfcN(N)/t*A8F 
GfcMCM)a(ifN(N)-BASEACARRY 
199  CONTINUE 

?oo  continue 
N^TARIaNSTRT 

IF  IN5TART.LE.0)  NST ARTa?uftl 
NSTARTai*(N5TART/2)4l 
00  300  Nai'NMKO 
NTEMPaNbTART/BABC 
RAN(N)aN$TART-NTEMP*»lA8fc 
300  NSTARTaNTcMR 
RETURN 
END 


Figure  B-2.  Fortran  listing  of  RANSET 
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START 


END 


Figure  B-3.  Logic  flowchart  for  URAND 
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Figure  B-4.  Logic  flow  chart  for  RANSET 
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i 

15 

GEN(IO)  -  An  array  containing  the  generator  X(=  5  )  in  multiple 
precision  representation 

NWRD  -  The  number  of  words  used  in  the  multiple  precision 

t 

representation  of  an  integer 

BASE  -  The  base  used  in  the  multiple  precision  representation 

MOD  -  The  maximum  value  of  the  highest  order  'digit'  in  the 

multiple  precision  representation 
FBASE  -  Floating  point  value  of  BASE 

FMOD  -  Floating  point  value  of  MOD 

RAN,  GEN,  NWRD,  and  NBASE  are  Fortran  integers;  FBASE  and  FMOD  are 
Fortran  real  quantities. 

An  alternative  method  (unfortunately,  not  machine  independent)  of  giving 
the  routine  a  starting  value  is  to  save  the  array  RAN  at  the  end  of  a  run  and  to 
restore  RAN  at  the  start  of  the  new  run  (just  after  the  RANSET  call). 

In  the  last  step  of  the  URAND  flow  the  objective  is  conversion 
of  the  multiple  precision  integer  random  number  R  to  a  floating  point 
random  number  X  between  0  and  1.  The  multiple  precision  integer 
produced  by  the  random  number  algorithm  is  represented  by  the  'digits' 
r^,  rg, . .  rn  (remember  that  r ^  is  the  lowest  order  digit.  Thus, 

R  =  r1  +  (BASE)-r2  + (BASE)2. r3 +....+ (BASE)N-1.rN  . 

Notice  that  we  have,  from  the  manner  in  which  N  and  MOD  were  established, 

r 

P  =  (BASE)**’’1.  MOD  . 
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The  uniform  random  number  desired  is  given  by  R/P.  Thus  we  have, 


X  = 


P  - - WIT - 

r  (BASE)”  A.  MOD 


(base)n'z-mod 


(base)n-S-  mod 


rN-l  rN 

*  •  •  •  55SB*  mod  *  KIOE 

=  M5D(rN  +  S?Bl(rN-l  +  ****  gffiSE(r2  +  BraF*rl)-***)) 
Starting  from  the  right  it  is  easy  to  compute  this  iteratively. 

B.  6  FIRST  100  RANDOM  NUMBERS  PRODUCED  BY  MIRAN 


For  checkout  purposes,  Table  B-l  lists  the  first  100  random  num¬ 
bers  produced  by  MIRAN  when  the  default  value  of  NSTART,  2001,  is  used 
as  the  starting  random  number. 
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REFERENCES  AND  ABSTRACTED 
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1.  Coveyou,  R.  R. ,  and  R.  D.  MacPherson,  ’’Fourier  Analysis 
of  Uniform  Random  Number  Generators,  ”  Journal  of  the  ACM, 

14  pp.  100-119,  1967. 

A  method  of  analysis  of  uniform  random  number  generators  is  de¬ 
veloped,  applicable  to  almost  all  practical  methods  of  generation. 

The  method  is  that  of  Fourier  analysis  of  the  output  sequences  of 
such  generators.  With  this  tool  it  is  possible  to  understand  that 
predict  relevant  statistical  properties  of  such  generators  and  com¬ 
pare  and  evaluate  such  methods.  The  results  of  many  such  analyses 
and  comparisons  are  given.  The  performance  of  these  methods 
as  implemented  on  differing  computers  is  also  studied  The  main 
practical  conclusions  of  the  study  are:  (a)  Such  a  priori  analysis 
and  prediction  of  statistical  behavior  of  uniform  random  number 
generators  is  feasible,  (b)  The  commonly  used  multiplicative 
congruence  method  of  generation  is  satisfactory  with  careful  choice 
of  the  multiplier  for  computers  with  an  adequate  35  bit)  word 
length,  (c)  Further  work  may  be  necessary  on  generators  to  be 
used  on  machines  of  shorter  word  length. 

2.  Kahn,  H. ,  Applications  of  Monte  Carlo,  Rand  Corp. ,  AEC-3259, 
USAEC,  April  1964. 

A  classic  publication  in  the  field  of  Monte  Carlo  methods  that  describes 
general  Monte  Carlo  methods,  random  number  generation  schemes 
and  variance  reduction  techniques.  The  volume  is  divided  in  two 
parts.  Part  I  describes  basic  techniques  with  random  numbers  (such 
as  fundamental  random  number  generation  techniques)  and  Part  n 
details  several  variance  reduction  schemes.  The  general  areas  of 
application  addressed  are  problems  in  radiation  transport. 

3.  MacLaren,  M.D. ,  G.  Marsaglia,  and  T.  A.  Bray,  ”A  Fast  Procedure 
for  Generating  Exponential  Random  Variables, "  Communications  of 
the  ACM,  7,  May  1964. 

A  very  fast  method  for  generating  exponential  random  variables  in  a 
digital  computer  is  outlines.  A  detailed  flow  diagram  and  required 
tables  are  provided. 
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4.  Marsaglia,  G,  and  T.  A.  Bray,  "A  Convenient  Method  for  Generating 
Normal  Variables,,  "  SIAM  Review,  6,  1964. 

A  ve.  */  fast  yet  small  Fortran  routine  for  generating  normal  random 
variables  in  terms  of  a  sequence  of  random  variables  uniform  rve r 
[  0,  l]  is  presented.  A  random  variable  X  is  generated  in  terms 
of  uniform  variables  IL ,  U«, . . . .  in  the  following  way:  86  percent 
of  the  time,  X=2(U i+uAu  f  -  1.  5),  11  percent  of  the  time,  X  =  1.  5 
(U.+Ug  -  1),  and  the  remaining  3  percent  uses  a  complicated  pro¬ 
cedure. 

% 

5.  Marsaglia,  G. ,  M.  D.  MacLaren,  and  T.  A.  Bray,  "A  Fast  Procedure 
For  Generating  Normal  Random  Variables, "  Communications  of  the 
ACM,  7,  1964. 

A  technique  for  generating  normally  distributed  random  numbers  is 
described.  It  is  faster  than  those  currently  in  general  use  and  is 
readily  applicable  to  both  binary  and  decimal  computers. 

6.  National  Bureau  of  Standards  Applied  Mathematics  Series  55,  June 
1964.  Handbook  of  Mathematical  Functions,  Numerical  Methods, 
pp.  849-953. 

This  section  of  the  handbook  reviews  various  methods  of  generating 
random  numbers  including  the  rejection  and  composition  methods. 

Also  presented  are  specific  techniques  for  various  discrete  and  con¬ 
tinuous  distributions  such  as  the  normal  and  exponential  distributions. 

7.  Spanier,  J. ,  and  E.  M.  Gelbard,  Monte  Carlo  Principles  and  Neu¬ 
tron  Transport  Problems,  Addision  Wesley  Ihiblishers,  1969.  "* 

This  is  one  of  the  more  recent  comprehensive  references  on  Monte 
Carlo  methods  as  applied  to  radiation  transport  problems.  Basic 
fundamentals  of  Monte  Carlo  are  first  reviewed.  Next  the  concepts 
of  discrete  and  continuous  random  walks  are  introduced  followed  by 
a  discussion  of  variance  reduction  techniques.  Finally,  advanced 
concepts  and  applications  to  radiation  transport  are  presented. 
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